Part 1

1.1 Blood of murderer

Высказывание прокурора: “Шанс, что у подсудимого была бы именно такая группа крови, если бы он был невиновен -- всего 1%; значит, с вероятностью 99% он виновен...“.

То есть он оценивал $P(has\_rare\_blood|not\_guilty)$, которую посчитал равной априорной вероятности выбрав произвольного человека, попасть в человека с редкой группой крови. Здесь под rare-подразумевается группой крови как у убийцы (она же редкая).

$$P(has\_rare\_blood|not\_guilty) = P(has\_rare\_blood) = \frac{\sum{people\_with\_rare\_blood\_group}}{\sum{all\_people}}=0.01$$

И здесь он уже не прав, потому что: $$P(A) = \sum{P(A|B_i)}$$

То есть: $$P(has\_rare\_blood) = P(has\_rare\_blood|not\_guilty) * P(not\_guilty) + P(has\_rare\_blood|guilty) * P(guilty) = P(has\_rare\_blood|not\_guilty) * P(not\_guilty) + 1 * P(guilty) = 0.01$$

Так как мы знаем, что $P(has\_rare\_blood|guilty)=1$, так как достовено известно, что у убийцы редкая группа крови. И следовательно:

$$P(has\_rare\_blood|not\_guilty) = \frac{P(has\_rare\_blood) - P(guilty)}{P(not\_guilty)}$$

И $P(has\_rare\_blood|not\_guilty) = P(has\_rare\_blood)$ лишь при предположении, что $\frac{P(guilty)}{1 - P(not\_guilty)} = 0.01$

При этом в суде должна интересовать другая вероятность, а именно: является ли данный человек убийцей (то есть виновен или нет), считая известным факт того, что у него редкая группа крови. Рассмотрим как меняется вероятность виновности человека, после того, как стало известно, что он является обладателем редкой группы крови, применив теорему Байеса

$$P(guilty|has\_rare\_blood) = \frac{P(has\_rare\_blood|guilty)}{P(has\_rare\_blood)} * P(guilty) = \frac{1}{0.01} * P(guilty)$$

1.2 Attorney's argument

Довод адвоката состоял в том, что он предположил, кровь могла равновероятно принадлежать любому человеку с такой редкой группой. И таким образом, он просто посчитал значение как отношение. (Хотя даже в такой постановке кажется странным отбрасывать довод, так как даже если нет никаких других улик, и выбор "виновного" равновероятен среди людей с группой такой крови, то 1e-4 гораздо лучше, чем ~1e-6).

В сущности же нужно рассматривать отношение улики к делу имено как отношение априорной вероятности вины к апостериорной (или относительной): $$\frac{P(guilty|has\_rare\_blood)}{P(guilty)} = \frac{1}{0.01} = 100$$

1.3 A terrible desease

Некоторые формулы: $$P(seak|positive\_test) = \frac{P(positive\_test|seak) * P(seak)}{P(positive\_test)}$$

Где: $$P(positive\_test|seak) = 1 - p_{false\_negative}$$ $$P(positive\_test|!seak) = p_{false\_positive}$$

$$P(seak) = p_{seak} = 0.01$$$$P(positive\_test) = P(positive\_test|seak) * P(seak) + P(positive\_test|!seak) * P(!seak) = (1 - p_{false\_negative}) * p_{seak} + p_{false\_positive} * (1 - p_{seak}) = 0.01 * (1 - p_{false\_negative}) + 0.99 * p_{false\_positive}$$

Подстваляя в формулу: $$P(seak|positive\_test) = 0.01 * \frac{1 - p_{false\_negative}}{0.01 * (1 - p_{false\_negative}) + 0.99 * p_{false\_positive}}$$

Аналогично: $$P(seak|!positive\_test) = \frac{P(!positive\_test|seak) * P(seak)}{P(!positive\_test)}$$

Где: $$P(!positive\_test|seak) = p_{false\_negative}$$ $$P(!positive\_test|!seak) = 1 - p_{false\_positive}$$

$$P(!positive\_test) = P(!positive\_test|seak) * P(seak) + P(!positive\_test|!seak) * P(!seak) = p_{false\_negative} * p_{seak} + (1 - p_{false\_positive}) * (1 - p_{seak}) = 0.01 * p_{false\_negative} + 0.99 * (1 - p_{false\_positive})$$

И финально: $$P(seak|!positive\_test) = 0.01 * \frac{p_{false\_negative}}{0.01 * p_{false\_negative} + 0.99 * (1 - p_{false\_positive})}$$

С точки зрения пользователя, мне бы хотелось уменьшить число ложно положительных срабатываний: меня бы не устраивала вероятность в 16% после получения положительного теста. С другой стороны, получение отрицательного результата теста дают мне шансы быть здоровым около 99,95%, что кажется вполне приемлимым результатом.

Part 2

2. Постройте графики целевых переменных.

Вы увидите, что число заболевших растёт очень быстро, на первый взгляд экспоненциально

2.a.

Используя линейную регрессию, обучите модель с экспоненциальным ростом числа заболевших: y ~ exp(линейная функция от x), где x — номер текущего дня.

2.b.

Найдите апостериорное распределение параметров этой модели для достаточно широкого априорного распределения. Требующееся для этого значение дисперсии шума в данных оцените, исходя из вашей же максимальной апостериорной модели (это фактически первый шаг эмпирического Байеса).

Вместо обучения целевой переменной $y$ - числа заболевших, будем обучать $z=ln (y)$, которая очевидно связана с целевой перменной. Кроме того, оценка на шум скорее выполняется именно для $z$ (о чем например, свидетельсвтуют графики функций, представленные выше.

В качестве априорного распределения возьмем также нормальное (сопряженно-априорное): $$p(w) = N(\mu_0, \sigma_0)$$ Из графиков выше понятно, что $w_0$ должно быть в районе $[0.2-0.3]$. Однако в задании было сказано выбрать достаточно широкое априорное распределение, поэтому можем взять, например $\mu_0=0, \sigma_0=1$:

Соответственно, прогнозируемая величина $z_i = ln (y_i)$. Будем исходить из предположения, что каждый y_i независим от предыдущих (хотя очевидно, что это не так). Тем не менее, оценим правдоподобие как: $$LH = p(data|w) = p(\overline{z}|w) = p(z_1,z_2,...z_n|w) = \prod_{i=1}^{n}p(z_i|w)$$ $$ln(p(data|w)) = \sum_{i=1}^{n}ln(p(z_i|w))$$

Будем исходить из предположения, что $z_i ~= f(x_i, w) + noise$ и при этом шум можно распределен нормально относительно предсказанной величины: $$z_i = f(x_i, w) + N(0, \Sigma) = N(f(x_i, w), \Sigma) = N(w^Tx, \Sigma)$$

Тогда апостериорная вероятность также является гауссианой: $$p(w|data) = \frac{p(data|w) p(w)}{p(data)} = C * p(data|w) p(w)$$ $$ln(p(w|data)) = + \sum_{i=1}^{n}ln(p(z_i|w)) = C' -\frac{nd}{2}ln(2\pi) - \frac{n}{2}ln(|\Sigma|) -\frac{1}{2}(X^Tw - \overline{z})^T\Sigma^{-1}(X^Tw - \overline{z}) - ln(2\pi) - \frac12ln(|\Sigma_0|) - \frac12(w-\mu_0)^T\Sigma_0^{-1}(w-\mu_0)=$$ $$= C'' - \frac12w^T(X\Sigma^{-1}X^T + \Sigma_0^{-1})w + (\overline{z}^T\Sigma^{-1}X^T + \mu_0^T\Sigma_0^{-1})w - \frac{n}{2}ln(|\Sigma|) - \frac12\overline{z}^T\Sigma^{-1}\overline{z}$$

Откуда собственно можно получить параметры апостериорного распределения: $$\Sigma_{posterior} = X\Sigma^{-1}X^T + \Sigma_0^{-1}$$ $$\mu_{posterior} = \Sigma_{posterior}(\overline{z}^T\Sigma^{-1}X^T + \mu_0^T\Sigma_0^{-1})$$

$$p(w|data) = N(\mu_{posterior}, \Sigma_{posterior})$$

Будем считать, что каждый сэмпл имеет одинаковый и независимый шум. В таком случае: $\Sigma = \sigma^2 I$

И в таком случае легко вычислить соответствующий максимальному правдоподобию шум: $$\frac{d(ln(p(w|data)))}{d\sigma^2} = \frac1{2\sigma^4}(X^Tw - \overline{z})^T(X^Tw - \overline{z}) - \frac{n}{2\sigma^2} = 0$$ $$\sigma^2_{max\_likelihood} = \frac1n(X^Tw - \overline{z})^T(X^Tw - \overline{z})$$

2.c.

Посэмплируйте много разных экспонент, постройте графики. Сколько, исходя из этих сэмплов, предсказывается случаев коронавируса в России к 1 мая 2020 года? к 1 июня? к 1 сентября? Постройте предсказательные распределения (можно эмпирически, исходя из данных сэмплирования).

Построим графики, аналогичные тем, что мы строили выше, только для большего числа дней:

Видно, что прогноз начинает плохо работать при отдалении от обучающих данных. Впрочем понятно, что экспоненциальный рост числа заболевших не мог продолжаться бесконечно долго.

3.

Кривая общего числа заболевших во время эпидемии в реальности имеет сигмоидальный вид: после начальной фазы экспоненциального роста неизбежно происходит насыщение. В качестве конкретной формы такой сигмоиды давайте возьмём форму функции распределения для гауссиана: $$\Phi(x) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{x}e^{-t^2/2}dt$$

Непонятно как обучать такую функцию. Но мы можем попровать обучить другую: $\Phi'(x)=e^{ax^2+bx+c}$ которая по смыслу является скоростью роста общего числа заболевших по дням, то есть фактически новое число заболевших в зависимости от дня.

Далее применим тот же подход и вместо обучения $y(x)$ - число заболевших в день $x$ будем обучать функцию $z(x)=ln(y)=w^T\overline{x}$, где $\overline{x} = \begin{pmatrix}x^2 \\ x \\ 1\end{pmatrix}$

Если построить график $z$, то становится видно, что лучше начать обучение с 12 марта (иначе 0 в начале слишком сбивают модель):

Вывод

Как видно из графиков, сигмоида гораздо лучше подходит для приближения результатов, ее она неплохо прогнозировала total_cases даже через месяц, после получения последней точки. Тем не менее, понятно, что она не может догадаться, когда будет пик, если данных ей не хватает. И уж тем более не может обнаружить "вторую волну" эпидемии

Попытка использовать не сопряженнно-априорные распределения и крах надежд:( (не читать)

Итак, предполагаем, что наша целевая переменная - распределена нормально относительно $e^{w^T * x}$, то есть $N(e^{w^T * x}, \sigma^2)$

Правдоподобие: $$p(t|w,\sigma^2,X) = \displaystyle\prod_{i=1}^{N} {N(e^{w^T * x_i}, \sigma^2)}$$ $$ln(p(t|w,\sigma^2,X)) = \displaystyle\sum_{i=1}^{N} {ln(N(e^{w^T * x_i}, \sigma^2))} = -\frac{N}{2}ln(2\pi\sigma^2) - \frac{1}{2\sigma^2}\displaystyle\sum_{i=1}^{N}(t_i - e^{w^T * x_i})^2$$

Найти $w$, удовлетворяющий условия максимума правдоподобия, можно попытаться, продифференцировав по $w$: $$\nabla_wln(p(t|w,\sigma^2,X)) = \frac{1}{2\sigma^2}\displaystyle\sum_{i=1}^{N}(t_i - e^{w^T * x_i}) * e^{w^T * x_i} * x_i = 0$$ Но решить такое аналитически мне сложно:)

Как выбрать "достаточно широкого априорного распределения" $p(w)$? Например, из графика понятно, что $w_1$ лежит в положительной части оси ($>0$), но при этом достаточно мал. Что же касается свободного члена ($w_0$), то относительно него сложно сделать какие то априорные предположения. Соответсвенно можно предложить использовать равномерное в некотором окне распределение на параметрах $w_0, w_1$. Или же какое-то убывающие на бесконечности, например нормальное относительно $w_0$ и "скошенное" нормальное (относительно $w_1$): $$w_0 = N(0, \sigma_{w_0}^2)$$ $$w_1 = \begin{cases} 0 & x\leq 0 \\ 2 * N(0, \sigma_{w_1}^2) & 0 < x \end{cases} $$

В таком случае апостериорное распределение параметров w: $$p(w|t,x) = \frac{p(t|w) * p(w)}{p(t)} = \frac{1}{const}p(t|w)*p(w) = \frac{2}{const} * sign(w_1)*N_{w_1}(0, \sigma^2_{w_1}) N_{w_0}(0, \sigma^2_{w_0})\displaystyle\prod_{i=1}^{N} {N(e^{w^T * x_i}, \sigma^2)}$$ Из нее можно выразить дисперсию шума: $$\nabla_{\sigma^2}ln(p(w|t,x)) = \nabla_{\sigma^2}\displaystyle\sum_{i=1}^{N}ln\frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{t_i-e^{w^T*x_i}}{2\sigma^2}} = \nabla_{\sigma^2}(-\frac{N}{2}ln (2\pi\sigma^2) - \frac{1}{2\sigma^2}\displaystyle\sum_{i=1}^{N}(t_i-e^{w^T*x_i})) = -\frac{N}{2\sigma^2} + \frac{1}{2\sigma^4}\displaystyle\sum_{i=1}^{N}(t_i-e^{w^T*x_i}) = 0$$ $$\sigma^2_{max\_likelyhood} = \frac{1}{N}\displaystyle\sum_{i=1}^{N}(t_i-e^{w^T*x_i})$$

Где в качестве $w$ используется $w_{max\_likelyhood}$: $$\nabla_wln(p(w|t,x)) = \nabla_w N_{w_1}(0, \sigma^2_{w_1}) N_{w_0}(0, \sigma^2_{w_0})\displaystyle\prod_{i=1}^{N} {N(e^{w^T * x_i}, \sigma^2)}$$ Значение константы в апостериорном распределении получается из интеграла: $$const = p(t) = \int_{-\inf}^{\inf}p(t|w)*p(w) dw$$